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Abstract. We study analytically and numerically a model describing front propagation 
of a KPP reaction in a fluid flow. The model consists of coupled one-dimensional reaction- 
diffusion equations with different drift coefficients. The main rigorous results give lower 
bounds for the speed of propagation that are linear in the drift coefficient, which agrees 
very well with the numerical observations. In addition, we find the optimal constant in a 
functional inequality of independent interest used in the proof. 



Many reaction processes in nature and engineering take place in a moving fluid. The 
creation of ozone in the atmosphere, the nuclear explosion of a supernova, pattern formation 
in morphogenesis, wild fires, and the gasoline transformation in an internal combustion 
engine are just a few examples. The effect of fluid flow on a reaction process can be 
profound, especially if the flow is strong. The problem has been studied for many years 
by engineers, physicists, and mathematicians alike. One of the most extensively studied 
mathematical models of the premixed reaction process is the reaction-diffusion equation 
and systems. The advection-reaction-diffusion equation is given by 



Here u(x) is the fluid flow, which we will assume is time-independent and passive. The 
coefficient A is the flow strength parameter, k is the diffusivity, M is the reaction strength, 
and f(T) is the reaction term. The function T(x,t) is normalized so that < T(x,t) < 1, 
and can stand for normalized temperature, the mass fraction of a reactant, or a share of 
a population with a certain feature, depending on the problem. Classical works of Fisher 
and Kolmogorov-Petrovski-Piskunov (KPP) ^21 Ej first considered equation (jl.lj) in the 
dimension d — 1 and with u — 0. These works modeled the propagation of an advantageous 
gene in the population, and established existence and, in a certain sense, stability of the 
traveling fronts. The extensions have occupied some of the best efforts of mathematicians 
since. Equation (jl.lj) in several spacial dimensions in the absence of advection (u = 0) is 
by now fairly well understood: the existence of traveling waves, stability, and asymptotic 
propagation properties have been extensively studied (see, e.g. EH EH 121] where many 
more references can be found). The effects of advection, however, are very important in 
many situations in combustion and more general chemical reactions fZT\ |3UJ |32] as well as 
in biology and ecology [221 1231 ■ 
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1. Introduction 




T t + Au ■ VT = kAT + Mf(T). 
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One important phenomena that is observed in many situations is the enhancement of the 
reaction rate by fluid motion. The physical reason for this observed speed-up is believed to 
be that fluid advection tends to increase the area available for reaction. Recently, there has 
been significant progress in the mathematical understanding of the reaction rate enhance- 
ment for several different classes of flows [H [21 SI O HOI EH! UHl UHl HZl [201 123 I2H| - 
The papers jSJ HH] provide excellent reviews of some of these results and further references. 

In this paper, our goal is to study reaction enhancement within a framework of a model 
given by a coupled system of one-dimensional advection-reaction-diffusion equations: 

(1.2) (Tj) t + A,{T 3 ) X - (Tj) xx = 7}(1 - Tj) + a(T^ - T 3 ) + a(T j+1 - Tj), 

j — 1, . . . , N. We assume that the model is periodic in j; that is, in the equation for Ti(x), we 
replace T (x) with T^(x), and in the equation for T^(x), we put T\(x) instead of T^+i(x). 
The initial data always satisfies < Tj(x, 0) < 1, and it is a simple consequence of the 
maximum principle that Tj(x,t) remains within these bounds for all times (see Section |2J). 
Our reaction term is the classical KPP |2Ij, and for simplicity we set k = M = 1 (these 
parameters are normalized by a simple rescaling of time and space). We assume that 
^2j =1 Aj = 0, which corresponds to the mean zero flow. Any nonzero mean is taken into 
account by switching to a moving system of coordinates; we are interested in a non trivial 
effect of front stretching rather than simple transfer with a constant speed. One can regard 
f)1.2|) as a model for reaction in a layered fluid, where different layers move with different 
speeds. The role of the parameter a > is to provide diffusive coupling between different 
layers, and its meaning will be further discussed in Section 0] Our main interest is in the 
case of the large coefficients Aj, where the influence of the drift is most pronounced. We 
consider the front-like initial data where Tj(x, 0) — > as x — > oo and Tj(x, 0) ^ 1 as 
x — > — oo exponentially fast (we will state the exact conditions in Section |2J). We adopt the 
following natural quantity (called the bulk burning rate in [5|) as the main measure of the 
reaction rate: 

N °9 

(1.3) V (t) = ^J2 / ( T M*)dx. 

J —oo 

One of the results that we prove is as follows. 

Theorem 1.1. Assume that the initial data Tj(x,0) are front like (|2.2j) . Then there exists 
a universal constant C such that for any r > 0, we have 



i r r N 

1.4) - V{t)dt>-'£ i \A j -. 1 -A j 

n ; =1 



1 1 

+ — + - + a + l 



ar 2 ar r 



There are several noteworthy properties the lower bound (jl.4j) possesses. Firstly, the 
bound is linear in the flow strength. Secondly, the estimate clearly shows the need to wait 
for a certain time r before the front propagation (and thus the reaction rate) stabilizes and 
the long-time lower bound is achieved. If the coupling constant a is large, we have the lower 
bound CAor x starting from times r ~ 1. If a. is small, the lower bound is CA, starting 
from times r w oT 1 . We interpret these results further in Sections 0] and El The numerical 
experiments discussed in Section El show very good qualitative agreement with ()1.4j) . 
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We remark that the result of Theorem 11.11 is reminiscent of the bounds proved in [§| for 
the two dimensional equation (jl.l|) with a shear flow u = (u(y),0). This is not surprising 
since our system (|1.2jl can be regarded as a discrete model of a shear flow. However, the 
proofs in the case of model (jl.2j) are conceptually much more transparent, the estimates are 
more precise, and the dependence on the averaging time is more explicit. 

The paper is organized as follows. In Section |21 we collect some background estimates 
which are needed in the proof of main results. In particular, we prove a sharp version 
of a functional inequality which was first discovered in Pj. The inequality is somewhat 
reminiscent of the indeterminacy principle of quantum mechanics. Although finding a sharp 
constant is not particularly important for the key results of this paper, we feel that the result 
is elegant (even if elementary) and may be of independent interest. In Section |H] we discuss 
a simpler two-layer case. In Section 0] we prove Theorem 11.11 and some related results for 
a slightly more general model. In Section El we describe the numerical simulations. 

2. The Auxiliary Tools and The Optimal Constant 

We make our model slightly more general by associating a width hj with each layer. By 
dimensionality of the discrete gradient terms providing coupling of the neighboring layers 
in (jl.2|) . the parameter a is made to depend on j, with ctj — n/hj : 

(2.1) (Tj) t + Aj(Tj) x - (Tj) xx = Tj{l - Tj) + a^T^ - Tj) + aj (T j+1 - Tj). 

The parameter k now plays the role of diffusivity between the layers. As opposed to the 
diffusivity in x direction, this diffusion coefficient cannot be set to one by a simple rescaling. 
Our mean zero flow condition now reads X)j=i fyj'A? = 0. Set H = J2f=i hj- 

We first prove a lemma on the preservation of spacial decay of the solutions of (|2.1|) that 
will allow us to manipulate the equations, in particular integrate by parts. 

Lemma 2.1. Assume that Tj(x,t), j = 1,...,N satisfy (J2.1j) and that the initial data 
Tj(x, 0) satisfy 

(2.2) T 3 (x,0)<C e- Xx , 1 - Tj{x, 0) < C e Xx , \(Tj{x, 0)) x \ < C e~ x ^ 
for all j and some A > 0. Assume that c\ is such that 

(2.3) C\ > m&Xj\Aj\ + A + A _1 (l + 2maXjOj). 
Then for all x, t > we have 

(2.4) Tfat) < C e~ xix - Clt) , 1 - Tjix^) < C e x{x+Clt \ \(Tj(x,t)) x \ < C e x{M+Clt) . 

Proof. Let us start by establishing the first bound in (J2.4|) ; the second bound can be proved 
by an identical argument. Set <f>{x,t) = e~ x( - x ~ Clt \ and Gj(x,t) = Co<f>(x,t) —Tj(x,t). A 
direct computation using ()2.3|) shows that 

(2.5) <p t + maXj\Aj\(j) x — <p xx — (1 + 2max :/ -(x ; -)0 > 

for any j. Note that Gj(x,0) > in view of (|2.2j) . On the contrary, assume that the first 
bound in (|2.4|) is not satisfied and that to > is the first time such that there exists j and 
Xq so that Gj(x ,t ) = 0. Combining ()2.5j) with (|2.1|) . we find that 

(2.6) ^ Gj ) t + Aj (G j ) x -(G j ) xx >0 
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for all x, t < t . Now, in ([2.6)1 . Gj(x,0) > for all x and Gj(x ,t ) = contradict the 
well-known maximum principle for parabolic PDE (see e.g. [llj). The proof of the second 
inequality in ()2.4j) is similar and is omitted. To prove the bound for (Tj) x , differentiate 
(12.1)1 : denoting Hj = (Tj) x , we have 

(".j)' + MHj) x - (Hj),-,- = U,(\ - ■>T J ) + 0,(7/, , - H s ) + a j+1 (H j+1 - II,). 

Given that < Tj{x) < 1 by a simple application of maximum principle, the bounds for 
(Tj) x are now obtained in the same manner as above. □ 

The bulk burning rate for the model ()2.1j) is now defined accordingly via 
(2.7) HO jjY.li, f (Tj)iUXl-r. 

7 = 1 J .. 



Assuming the initial data satisfies ([2.2)1 . the indefinite integrals are well-defined (integrating 
over [—B,B] and taking B to infinity). As a consequence of (|2.1)1 . Lemma 12.11 and the 
condition fyjA? = 0, we also have 



(2.8) V(t) = ^ / - T s {x)) dx. 

i= l -co 

We continue by computing: 

N OO OO 

^'(t) = ^X)/i i / r i (x)(l-T J -(x))dx=i^/i J - f(T j ) t (x)(l-2T j (x))dx 

J — oo J — oo 

^ — ^ — oo 

(2.9) > 1 £ ^ | (2((T,),) 2 + ^(T,(x) - Tj_x(x)) 2 ^ rfx - V(f). 

J' - — oo 

In the second step above, we substitute (Tj) t (x) from ()2.1j) . integrate by parts, and rearrange 
the terms. These manipulations are justified by Lemma 12.11 In the last step we use the 
expression ()2.8)) and the fact that < Tj(x) < 1. 

As a warm up before our main results, we prove the following estimate. 

Theorem 2.2. For any choice of Aj such that J^. hjAj = 0, we have 

TV 2 ( TV 2 



(2.10) V(tY>- + e- 2t [V(Q) i( . 

A key element in the proof is a general functional inequality. 
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Proposition 2.3. For any function T(x) € C 1 (M), < T(x) < 1, such that T(x) — > 1 as 
x — > — oo and T(x) — > as x — > oo, we nave 

oo oo 

(2.11) I (T x .) 2 rfx I T(l-T)dx> Q\ 

— oo — oo 

T/ie constant (|) 2 zs sharp. 

Proof. We note that the inequality (J2.ll)) has been proved in jS] with a weaker universal 
constant on the right hand side. Using the Cauchy-Schwartz inequality we find that 

oo oo / oo \ 2 

J {T x fdx j T(l-T)dx > I J T x (T(l-T)) l/2 dx 

— oo — oo \— oo 

Changing the variable and taking into account the asymptotic behavior of T(x), we have 

oo 1 

J T x (T(x)(l - T{x))) l ' 2 dx = J (T(l - T)f/ 2 dT. 

-oo 

The latter integral is computed explicitly by setting T = cos#, and evaluates to |. This 
proves (J2.ll)) . To show that this constant is sharp, note that if T x = (T(l — T)) 1//2 , then 



2 

7TN 2 



J (T x fdx j T(l-T)dx= j T X {T{1 — T)) l l 2 dx J = ( 

— oo — oo \-oo / 

We can solve for T(x) explicitly and use the solution to construct an explicit function 
satisfying the necessary asymptotic behavior for which the equality in (|2.11jl is achieved. 
In particular, 

[ 1 if x < -f 

T(x) = I l=fp if - f < x < f 
{ if f < x. 

is one such function. □ 
Proof of Theorem \2. 6 A From ([2.9)1 . we see that 



h*)+v(*)>|e^ / m) x fdx 
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From (|2.2|) and Proposition 12.31 we derive 



2 N U 

hj_ 

32H z - 00 



> 



3=1 J 7}(a;)(l — dx 

— oo 

x ^ 3 J 3 3 j 32V (t) 



,, / :V \ 1 2 



On the second step we use Cauchy-Schwartz inequality. Therefore, we have 

The inequality (|2.10j) follows from the standard application of the Gronwall lemma. □ 

Theorem 12 . 21 gives a universal lower bound on the reaction rate, independent of the flow. 
One consequence is that in the framework of our model, a mean zero flow cannot quench the 
reaction. However, we are interested in the bounds that show enhancement of the reaction. 
We will prove such bounds in the next two sections. 

3. The Two-Layer Model 

For the sake of clarity, we start with a simplified version of our general model (j2.1j) where 
only two layers of unit width are involved. Namely, we look at a system 

(3.1) (Tx) t + A(Tx) t - {Ti) xx = Tx(l - T x ) + k(T 2 - T x ) 

(T 2 ) t - A{T 2 ) t - (T 2 ) xx = T 2 (l - T 2 ) + k{T x - T 2 ) 

with the initial data T 12 (a;, 0) satisfying (|2.2|) . We set h\ — h 2 — 1 in the definition ()2.7|) of 
the bulk burning rate. Throughout the rest of the paper we denote by C different universal 
constants which enter the estimates. 

Theorem 3.1. There exists a universal positive constant C such that for any r > 0, we 
have 



(3.2) -j- / V{t)dt > CA 

o 



1 1 



KT 2 KT T 



K + l 



Remark. Note that Theorem 13. II is a particular case of Theorem II .11 for N = 2. 



Proof. We first integrate (j3.1|) over all real axis. Given the asymptotic behavior of T± and 
T 2 , ()2.2|) and Lemma \2. 11 we obtain 

(3.3) k(Ti)tdx — A — J M Ti(l — Tx)dx + « / R (T 2 - T x )dx 

(3.4) UT2)tdx + A = / R T 2 (l - T 2 )dx + k J^T, - T 2 )dx. 
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Subtracting (|3.4|) from ()3.3|) and integrating the result over t £ [t±, t 2 ], we find that 

{T x {r 2 )-T 2 {r 2 ))dx- [ (T^n) - T 2 (n))dx 



(3.5) 




+ / / (T 2 (l - T 2 ) - T\(l - T x ))dxdt 




+ 2k / {Ti- T 2 )dxdt = 2A(r 2 - n). 



We now claim the following Lemma. 
Lemma 3.2. For any t > 0, we have 



(3.6) 



/ \T x (x,t) -T 2 {x,t)\dx < 3 
Jr 



(T x (x, t)(l — Ti(x, t)) + T 2 (x, t)(l — T 2 (x, t))) dx 



(Ti(x,t) — T 2 (x,t)) 2 dx+ 



<-(V(t)+V'(t)) + 3V{t). 

K 



Proof. We begin by proving the first inequality, which is valid for any functions T x {x) 
To(x) E L 1 (M) fl L 2 (R) : 



(3.7) 



\T X -T 2 \dx < 3 



/ (T x - T 2 fdx + / (T t (l - TO + T 2 (l - T 2 ))dx 



Consider two different cases. Denote S the set of all x such that \Ti(x)—T 2 (x)\ > 1/3. Then 
\Ti(x) — T 2 (x)\ < 3\Ti(x) — T 2 (x)| 2 and ()3.7j) holds with C = 3 if we restrict the integration 
to set S. Next, assume that \Ti(x) — T 2 (x)\ < 1/3. Note that if at least one of Ti(x), T 2 (x) 
belongs to the interval [1/3,2/3], then T x {x){\ - T x {x)) +T 2 (x)(l -T 2 {x)) > 2/9. Consider 
the alternative where both T x (x) and T 2 (x) lie in either [0, 1/3) or (2/3, 1]. Note that the 
absolute value of the derivative of function x(l — x) satisfies |(1 — 2x)| > 1/3 in these 
intervals. Then by the mean value theorem, 

T 1 (x)(l-T 1 (x))+T 2 (x)(l-T 2 (x)) > |T 1 (x)(l-T 1 (x))-T 2 (x)(l-T 2 (x))| > i|Ti(x)-T 2 (a:)|. 

Therefore, if x £ R \ S, then 

|Tx(x) - T 2 (x)\ < 3 (7i0r)(l - T x (x)) + T 2 (x)(l - T 2 (x))) . 

Combining the two cases, we obtain ()3.7j) . 

Now the second inequality in (J3.6|) follows from the definition (J2.7j) of bulk burning rate, 
V(t), and the estimate (j2.9|) (recall that we set h X)2 = 1). □ 

Given Lemma [3.21 (J3.5|) implies that 
(3.8) -[V'(r 2 ) + V'in)] + (- + 9)[V(r 2 ) + V(n)} + (6k + 7) I V(t)dt > 2A(r 2 - n ). 

K K J 



Tl 
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Let us now apply the following averaging over t G [0, r] to the difference of equations 
Q and (jH3J): 



(3.9) 

where 
(3.10) 




dtdkds 



1 ' G(-,t--)dt 



l-sl-k 



G(h,0 



K^-iei) 2 -(|-iei) 2 if o<iei<| 



if § < iei < h. 



After the first integration, we obtain (|3.5|) and thus (|3.8|) with T\ = r — k, r 2 = r + /c. 
Integrating twice more leads to the following estimate: 



(3.11) 
where 



T T 

^ j V(t)dt +0 + 9^ y ^(t)/i(t, r)d* + 


T 

/t t At 3 

^M-,t--)rft> — 



h(t,r) 



t if o<t<i 

I if § < t < f 
r - t ii \ <t<r. 



Observe that G(^,t — |) < ^ and that h(t,r) < j. Hence (|3.11|) can be expressed as 



follows: 



T 

~ y v^(t)^ > ca 



111 

+ — + - + K + 1 



KT 2 KT T 



This coincides with (13 .2|) . thus proving the theorem. It is also easy to estimate that the 
constant C can be taken equal to 1/48. □ 

4. The Multi-Layer Model 

Here we prove a more general version of Theorem 11.11 Let us introduce the notation 

Ri = 1 + maxjhj 1 , R 2 = 1 + maxj^hjhj}^). 

Theorem 4.1. Assume that the initial data Tj(x,0) for (j2.1|) are front like ()2.2|) . Then 
there exists a universal constant C such that for any r > 0, we have 

N r . i-l 



(4.1) 



1 { vm>-^Y.\ A i- 



HR 2 



1 - A j\hj 



1 + — + — + kR 1 + 1 



KT KT T 



Remark. Theorem II .11 is a particular case of Theorem 14. II when hj = 1 (with a substituted 
by «)■ 
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Proof. Consider two neighboring layers, j and j — 1. Assume without loss of generality that 
Aj > Aj-i. Subtracting the equation for Tj from the equation for Tj_i and integrating, we 
find: 

(4.2) f ((Tj_ x ) t - (Tj) t ) dx + A, - Aj_i = [ (T^(l - T^) - Tj(l - Tj)) dx + 
Jm. Jr 

+ 7 / ( T i ~ T i-i) dx+ — (Tj - T j+1 ) dx + / (7)_ 2 ~ dx. 

hj-ij Jr n j Jr %-i Jr 



hj 



Multiplying (|4.2|) by hj, and summing over j — 1, . . . , N, we obtain: 

N „ N 

J^i-iphj / (^-i - Tj)tdx + hjlAj-t - Aj\ = 

N N , 

-J2(-lphj / Tj(l-Tj)dx + J^i-lphj^ / Tj^l-Tj^dx-^- 
j=l Jr j=1 Jr %-i 



(4.3) 



+ Ti-lpn \(l + p-) [ (Tj - Tj_ x )dx + / (Tj - T J+1 )dx 
~[ LV n j-i/ Jr Jr 



(Tj-2 - Tj_ t )dx 



h j-i 
hj_ 
hj-i 



where o~j = if Aj > Aj_\ and o~j = 1 otherwise. Denote the last sum in f)4.3|) by I\. Usinj 
Lemma f3. 21 and (|2.9jl . we can estimate: 



(4.4) 



III < Ck V / (Tj - T, x fdx + / T,(l - T,) f 1 + max (p- 
Ur Jr J V j \hj-i 



< CH (R 2 V'(t) + KV^R-yRi) . 



As in the proof of Theorem 13 ,1[ we average in i G [0, r] according to ()3.9|) and apply 
Lemma f3. 21 to estimate the expressions arising from the time derivatives in ()4.3)1 in terms of 
V and V. Taking into account (|4.4|) . after the first integration we obtain, with n = | — k 
and r 2 = | + fc, 



X>1)<^ / (2}_i(t 2 ) - Tj(r 2 ))dx - [ (Tj^(n) - Tj(r x ))dx 
j =1 Ur Jr 

(t 2 - n) ^ ^-l^-i - A- 1 < c™ 2 I V(r 2 ) - l/(n) + + l) y V(t) dt J . 
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Applying Lemma 13.21 once again to the expression on the left hand side, and integrating 
twice more, we obtain 



T -Y,hj\Aj^-Aj\<CR 2 



N 



3=1 



T 

Mi + 1) J V(t)G(^,t- T -)dt+ 



T T 

[k" 1 + J V(t)h(t,T)dt + k- 1 J V(t) dt 





Taking into account the bounds on functions G and h, we arrive at 

(4.5) - [ V(t) dt > V hjlAj-i -A j \(l + KR 1 + ^- + — + —) . 

T J HR 2 ^ V KT KT T J 

3=1 

This completes the proof of Theorem 14.11 In the situation where all hj = 1, and hence 
Ri = R2 = 2, we obtain Theorem 11.11 exactly. In general, we see that small hj (which 
translate to the large coupling constants ctj), or strongly non-uniform size distributions of 
the neighboring hj, weaken the lower bound. 

□ 



5. Numerical Simulations 

We carried out numerical simulations for the two layer system (|3.1|) . Our main interest 
was in finding the reaction propagation rate in the established regime, depending on A and 
k. Our results are obtained using the Crank-Nicholson implicit scheme. The implementation 
features two noteworthy modifications. Firstly, the transient wave that occurs between the 
initial condition and the eventual steady wavefront shape was ignored when calculating 
the front speed. This was achieved by using a " t transient " parameter, and then ignoring the 
wave speed for times to transient- Secondly, we employed a computational trick to expedite 
simulation. In a naive implementation, we would apply the Crank-Nicholson scheme to the 
entire field of points that the wave travels through. Since the matrices involved in the 
computation are n x n where n is the number of points being considered, considering the 
entire field of points that the wave passes through can be computationally expensive. We 
avoided this by only applying Crank-Nicholson to points in the immediate neighborhood 
of the front itself. Then after every iteration we adjusted our calculation frame to stay 
centered on the wavefront. 

We found clear linear dependence on the flow strength parameter A. The Figured] shows 
this linear dependence for different values of the coupling constant k. The constant k varies 
from 0.1 to 1 with 0.1 step. The slope of the graph is monotone decreasing in k, so that the 
steeper slopes on Figure ^ correspond to smaller values of k. The intuitive explanation for 
this effect is that stronger diffusion between layers mollifies the front stretching produced 
by the flow, thus reducing the reaction zone. The graph on Figure El shows the dependence 
of slope on k. The shape of the graph is in good qualitative agreement with the bound of 
Theorem 13.11 We note that the results of our simulations are similar to the results of 
where combustion in a shear flow was studied in a full PDE setting by means of a more 
sophisticated numerical scheme. 
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Wavefront speed vs advection 




Figure 1 . The bulk burning rate versus flow strength A for different values 
of the coupling k 

Slopes vs Coupling 
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Figure 2. The slopes corresponding to different values of k 
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